**************************************
*Replication code for 
*COVID-19 vaccination uptake in remote areas – Evidence from a panel survey in Bangladesh
*
*Published in PLOS One
*
*written by Lukas Rudolph
*
*2024-06-06
*
*contact: lukas.rudolph@uni-konstanz.de
*
*Data used: BD_df_covid_replication.dta 
*Data output generated: none
*Tables generated: Main Article Fig 2-6; Supplementary Information Fig 2-4
*Figutes generated: Main Article Table 1-3; Supplementary Information Tables 1-2
**************************************


********************************************************
*prepare stata
********************************************************

*set working directory to path where data is
cd "" // fill in path

*install packages if necessary: uncomment, find and install
/*
findit grstyle // by Ben Jann
findit estout // by Ben Jann
findit coefplot // by Ben Jann
findit tabplot // by Nicholas Cox
findit stripplot // by Nichloas Cox
*findit estpost // by Ben Jann
*/

*set grstyle for nice figure layout
grstyle init
grstyle set imesh
grstyle set legend, nobox

********************************************************
*load data
********************************************************

use BD_df_covid_replication.dta, clear // dataset reduced to variables necessary to replicate the analyses in this manuscript, and to sample of household heads only


******************************************************************************************************

* TABLES MAIN PAPER

******************************************************************************************************

**********************************
*Table 1
**********************************

eststo clear
estpost summarize agecat_tab? incomesource2_tab? Q3_demo_gender_tab? Q_demo_education2_tab? Q_demo_religion_enc_tab? below15_tab? above65_tab? w5_vacc  threat_family_health health_satisfaction any_affectedness end_expected  wave_expected susceptibility1 susceptibility2 village_vacc village_affected capital_duration_hrs notavailable ///
 ///
if V3_wave == 5 , detail
esttab, cells("mean(fmt(4)) p50(fmt(2)) p75(fmt(2)) p25(fmt(2)) count(fmt(0))") nomtitle nonumber label replace

*table 1 reformatted into form as presented in paper ; IQR = p75-p25

**********************************
*Table 2
**********************************

*set variables which are core components of HBM
global health_belief_model = "threat_family_health health_satisfaction i.any_affectedness end_expected  wave_expected susceptibility1 susceptibility2 village_vacc village_affected "

eststo clear

eststo: reg w5_vacc $health_belief_model  capital_duration_hrs notavailable i.V4_location /// 
 , cluster(V4_location)    
   vif // 3.60, non over 5 but capital duration hrs (34.99)

eststo: reg w5_fullvacc  $health_belief_model  capital_duration_hrs notavailable i.V4_location /// 
 , cluster(V4_location)   
   vif // 3.60, non over 5 but capital duration hrs (34.99)
 
 eststo: reg w5_vacc  $health_belief_model  i.V4_location /// 
 , cluster(V4_location)   
    vif // 1.77, non over 5  

eststo: reg w5_fullvacc $health_belief_model  i.V4_location /// 
 , cluster(V4_location)  
    vif // 1.77, non over 5  

	
  esttab  ///  using "$tab\tab2.tex" /// uncomment "using ..." in to save table
, label refcat (threat_family_health "\textbf{Health motivation}" end_expected "\textbf{Severity}" susceptibility1 "\textbf{Susceptibility}" village_vacc "\textbf{Social cues}"  capital_duration_hrs "\textbf{Availability}"  , nolabel) ///
indicate("Location fixed effects" = *V4_location ) ///
nobaselevels ///
b(2) se(2) star(+ 0.10 * 0.05 ** 0.01 ** 0.001) stat( N r2, fmt(0 2)) replace  nonote

 
 
**********************************
*Table 3
**********************************

global demographics = " i.agecat i.incomesource2 i.Q3_demo_gender i.Q_demo_education2 i.Q_demo_religion_enc i.below15 i.above65 " 
global fes = "i.V4_location"

eststo clear

eststo: reg w3_vacc    $demographics $fes ///
 , cluster(V4_location)
 vif // 1.75, non over 5
 eststo: reg w4_vacc    $demographics $fes ///
 , cluster(V4_location)
  vif // 1.72, non over 5

 eststo: reg w5_vacc    $demographics $fes ///
 , cluster(V4_location)
   vif // 1.70, non over 5

 eststo: reg w6_vacc    $demographics $fes ///
 , cluster(V4_location)
   vif // 1.72, non over 5

 eststo: reg w10_booster $demographics $fes ///
 , cluster(V4_location)
  vif // 1.74, non over 5

 
 esttab    ///   using "$tab\tab3.tex" /// uncomment "using ..." in to save table
, label refcat (notavailable "Availaibility" 1.agecat "\textbf{Age}" 1.incomesource2 "\textbf{Income source}" 1.Q2_demo_marital2 "\textbf{Marital}" 0.Q3_demo_gender "\textbf{Gender}" /// 
0.Q_demo_education2 "\textbf{Education}" 1.Q_demo_religion_enc "\textbf{Religion}" 0.below15 "\textbf{Children}" 0.above65 "\textbf{Elderly}" 0.affectedness "\textbf{Affectedness}" health_satisfaction "\textbf{Health satisfaction}" Q_trust_nat_gvmt_competence "\textbf{Gov. comp.}" 0.Q_attitudes_party_support "\textbf{Party ID}" 0.Q118_infec_infected "\textbf{Covid}" 1.vacc_village "\textbf{Village vacc.}" 1.V6_district_enc "\textbf{District}" , nolabel) ///
indicate("Location fixed effects" = *V4_location ) ///
b(2) se(2) star(+ 0.10 * 0.05 ** 0.01 *** 0.001) stat( N r2, fmt(0 2)) replace  nonote
 

******************************************************************************************************

* FIGURES MAIN PAPER

******************************************************************************************************

**********************************
*Figure 2
**********************************

* sample size
tab personal_infected V3_wave 
tab Q120_infec_infected_village_cnt V3_wave 

*figure
tabplot Q120_infec_infected_village_cnt V3_wave , percent(V3_wave) showval(format(%2.1f)) ytitle("How many in village affected?") ylabel(6 "None" 5 "Few" 4 "<50%" 3 "~50%" 2 ">50%" 1 "All")   subtitle("") name(one, replace)

tabplot personal_infected V3_wave , percent(V3_wave) showval(format(%2.1f)) ytitle("Personally affected?") ylabel(4 "No" 3 "Symptoms" 2 "Pos. test" 1 "Hospitalized")   subtitle("") name(two, replace)

graph combine two one, note(Percent of respondents by wave. N: 1'362-1'594,  size(*0.8))

*graph export "Fig2.eps", replace // uncomment to save figure



**********************************
*Figure 3
**********************************

*sample size 
tab Q121_vacc_vaccinated V3_wave

*figure
tabplot Q121_vacc_vaccinated V3_wave, percent(V3_wave) showval(format(%2.1f)) subtitle("") name(one, replace) ylabel(5 "No" 4 "Once" 3 "Twice" 2 "Thrice" 1 "4 times") xlabel(1 "Sept 21" 2 "Oct 21" 3 "Dec 21" 4 "Jan 22" 5 "Oct 22")

tabplot Q130_vacc_vaccinated_village V3_wave, percent(V3_wave)   yla(6 "None" 5 "Few" 4 "<50%" 3 "~50%" 2 ">50%" 1 "All" , ang(h)) showval(format(%2.1f)) subtitle("") name(two, replace) xlabel(1 "Sept 21" 2 "Oct 21" 3 "Dec 21" 4 "Jan 22" 5 "Oct 22")

graph combine  one two, note(Percent of respondents by wave. N: 1'327 -1'523,  size(*0.8))

*graph export "Fig3.eps", replace // uncomment to save figure



**********************************
*Figure 4
**********************************

*sample size 
tab Q125_vacc_AcceptVacc V3_wave, col

*figure
twoway histogram Q125_vacc_AcceptVacc if V3_wave==3 | V3_wave==5, by(V3_wave, note("") xrescale row(1)) horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(0 "No" 1 "Yes" 2 "Depends", ang(h)) discrete  xlabel(0(20)100)

*graph export "Fig4.eps", replace // uncomment to save figure



**********************************
*Figure 5
**********************************

*sample size
tab Q123_vacc_vaccinated_reason1 V3_wave if V3_wave==3 | V3_wave==5, nol col

*figure
twoway histogram Q123_vacc_vaccinated_reason1 if V3_wave==3 | V3_wave==5, by(V3_wave, note("") xrescale row(1)) horizontal width(1) percent  subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Protect me" 2 "Protect family" 3 "Protect vulnerable" 4 "Protect community" 5 "Protect society" 6 "Mandatory" 7 "Other", ang(h)) discrete

*graph export "Fig5.eps", replace // uncomment to save figure


**********************************
*Figure 6
**********************************

*prepare
label define incomesource2 1 "agri." 3 "busin." 4 "empl." 5 "dayl." 6 "gvmnt" 11 "oth."
label values incomesource2 incomesource2


*define macros
global demographics = " i.agecat i.incomesource2 i.Q3_demo_gender i.Q_demo_education2 i.Q_demo_religion_enc i.below15 i.above65 "  
global fes = "i.V4_location"

*to correct spacing/labeling
replace V3_wave = 8 if V3_wave==10

*xtset
xtset V1_respondent_number V3_wave

*regress and plot
eststo a: xtreg Q121_vacc_vaccinated i.V3_wave##i.Q_demo_education2 i.V3_wave##i.agecat i.V3_wave##i.incomesource2 i.V3_wave##i.Q_demo_religion_enc i.V3_wave##i.below15 i.V3_wave##i.above65  i.V4_location , cluster(V1_respondent_number) 
margins , at(V3_wave = (3 4 5 6 8) Q_demo_education2 = (0 1 2 3) ) post
marginsplot , name(one, replace) ylabel(0(.5)3) xlabel(3 "09/21" 4 "10/21" 5 "12/21" 6 "01/22" 8 "10/22") ytitle("") title("") xtitle("")

eststo b: xtreg Q121_vacc_vaccinated i.V3_wave##i.Q_demo_education2 i.V3_wave##i.agecat i.V3_wave##i.incomesource2 i.V3_wave##i.Q_demo_religion_enc i.V3_wave##i.below15 i.V3_wave##i.above65 i.V4_location  , cluster(V1_respondent_number) 
margins , at(V3_wave = (3 4 5 6 8) agecat = ( 1 2 3 4) ) post
marginsplot , name(two, replace)  ylabel(0(.5)3) xlabel(3 "09/21" 4 "10/21" 5 "12/21" 6 "01/22" 8 "10/22") ytitle("") title("") xtitle("")

eststo c: xtreg Q121_vacc_vaccinated i.V3_wave##i.Q_demo_education2 i.V3_wave##i.agecat i.V3_wave##i.incomesource2 i.V3_wave##i.Q_demo_religion_enc i.V3_wave##i.below15 i.V3_wave##i.above65  i.V4_location , cluster(V1_respondent_number) 
margins , at(V3_wave = (3 4 5 6 8) incomesource2 = ( 1 3 4 5 6 11) ) post
marginsplot , name(three, replace) legend(rows(2))  ylabel(0(.5)3) xlabel(3 "09/21" 4 "10/21" 5 "12/21" 6 "01/22" 8 "10/22") ytitle("") title("") xtitle("")

eststo d: xtreg Q121_vacc_vaccinated i.V3_wave##i.Q_demo_education2 i.V3_wave##i.agecat i.V3_wave##i.incomesource2 i.V3_wave##i.Q_demo_religion_enc i.V3_wave##i.below15 i.V3_wave##i.above65  i.V4_location , cluster(V1_respondent_number) 
margins , at(V3_wave = (3 4 5 6 8) Q_demo_religion_enc = (1 2) ) post 
marginsplot , name(four, replace) legend(rows(2))  ylabel(0(.5)3) xlabel(3 "09/21" 4 "10/21" 5 "12/21" 6 "01/22" 8 "10/22") ytitle("") title("") xtitle("")

graph combine one two three four, row(2) xsize(10) ysize(14)


*graph export "Fig6.eps", replace // uncomment to save figure


 

******************************************************************************************************

* FIGURES SUPPLEMENTARY INFORMATION

******************************************************************************************************


**********************************
*SI Fig 2
**********************************

*sample size
tab1 Q134a_threat_family_health Q134c_threat_society_health Q134a_threat_family_economic  Q134c_threat_society_economic

*xtset
xtset V1_respondent_number V3_wave 

*figure
stripplot  Q134a_threat_family_health Q134c_threat_society_health Q134a_threat_family_economic  Q134c_threat_society_economic ///
, by(V3_wave, rows(1) note("")) cumul cumprob box refline vertical centre ///
xlabel(1 "Fam. health" 2 "Soc. health" 3 "Fam. econ." 4 "Soc. econ.", ang(45)) ytitle("Is Covid a serious threat to...?") ylabel(1 "Str. agree" 2 3 4 5 "Str. disagree") 

*graph export "SIFig2.eps", replace  // uncomment to save figure


**********************************
*SI Fig 3
**********************************

twoway histogram Q140a_threat_family_tb if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(one, replace) xlabel(0(20)100) ytitle("Family / Tb.") 
twoway histogram Q140c_threat_society_tb if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(two, replace) xlabel(0(20)100) ytitle("Village / Tb.")
twoway histogram Q140b_threat_village_tb if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(three, replace) xlabel(0(20)100) ytitle("Country / Tb.")
twoway histogram Q141a_threat_family_flood if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(four, replace) xlabel(0(20)100) ytitle("Family / flood")
twoway histogram Q141b_threat_village_flood if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(five, replace) xlabel(0(20)100) ytitle("Village / flood")
twoway histogram Q141c_threat_society_flood if V3_wave==5, horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Smaller" 2 "Similar" 3 "Larger", ang(h)) discrete name(six, replace) xlabel(0(20)100) ytitle("Country / flood")
graph combine one two three four five six, col(1) ysize(12)

*graph export "SIFig3.eps", replace  // uncomment to save figure


**********************************
*SI Fig 4
**********************************

twoway histogram Q131_vacc_MandatoryVacc_general if V3_wave==3 | V3_wave==5, by(V3_wave, note("") xrescale row(1)) horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Voluntary" 2 "Mandatory", ang(h)) discrete name(one, replace) xlabel(0(20)100) ytitle("General population")
twoway histogram Q131_vacc_MandatoryVacc_byJob if V3_wave==3 | V3_wave==5, by(V3_wave, note("") xrescale row(1)) horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Voluntary" 2 "Mandatory", ang(h)) discrete name(two, replace) xlabel(0(20)100) ytitle("Certain occupations")
twoway histogram Q131_vacc_MandatoryVacc_byGroup if V3_wave==3 | V3_wave==5, by(V3_wave, note("") xrescale row(1)) horizontal width(1) percent subtitle(, fcolor(ltblue*0.2)) bfcolor(ltblue) blcolor(blue*0.5) yla(1 "Voluntary" 2 "Mandatory", ang(h)) discrete name(three, replace) xlabel(0(20)100) ytitle("Vulnerable populations")
graph combine one two three, col(1) l2title("Should vaccination be mandatory for... ")

*graph export "SIFig4.eps", replace  // uncomment to save figure
 

 
******************************************************************************************************

* TABLES SUPPLEMENTARY INFORMATION

******************************************************************************************************


**********************************
*SI Table 1
**********************************

eststo clear

eststo: logit w5_vacc  $health_belief_model capital_duration_hrs notavailable i.V4_location /// 
 , cluster(V4_location)    
 
eststo: logit w5_fullvacc  $health_belief_model capital_duration_hrs notavailable i.V4_location  /// 
 , cluster(V4_location)   
 
 
 eststo: logit w5_vacc  $health_belief_model i.V4_location /// 
 , cluster(V4_location)   
 
eststo: logit w5_fullvacc  $health_belief_model i.V4_location /// 
 , cluster(V4_location)  

eststo: logit w5_vacc  $health_belief_model capital_duration_hrs notavailable   /// 
 , cluster(V4_location)    
 
eststo: logit w5_fullvacc  $health_belief_model capital_duration_hrs notavailable  /// 
 , cluster(V4_location)   
 
 eststo: logit w5_vacc  $health_belief_model  /// 
 , cluster(V4_location)   
 
eststo: logit w5_fullvacc  $health_belief_model  /// 
 , cluster(V4_location)  
 
  
  esttab   ///   using "SITab1.tex" /// uncomment to save table
, label refcat (threat_family_health "\textbf{Health motivation}" end_expected "\textbf{Severity}" susceptibility1 "\textbf{Susceptibility}" village_vacc "\textbf{Social cues}"  capital_duration_hrs "\textbf{Availability}"  , nolabel) ///
indicate("Location fixed effects" = *V4_location ) ///
nobaselevels ///
b(2) se(2) star(+ 0.10 * 0.05 ** 0.01 ** 0.001) stat( N r2_p, fmt(0 2)) replace  nonote

 
**********************************
*SI Table 2
**********************************

*set variables for regression
global demographics = " i.agecat i.incomesource2 i.Q3_demo_gender i.Q_demo_education2 i.Q_demo_religion_enc i.below15 i.above65 " 
global fes = "i.V4_location"

eststo clear

eststo: logit w3_vacc    $demographics  $fes ///
 , cluster(V4_location)

 eststo: logit w4_vacc    $demographics  $fes ///
 , cluster(V4_location)

 eststo: logit w5_vacc    $demographics  $fes ///
 , cluster(V4_location)

 eststo: logit w6_vacc    $demographics  $fes ///
 , cluster(V4_location)

 eststo: logit w10_booster $demographics  $fes ///
 , cluster(V4_location)
 
 eststo: logit w3_vacc    $demographics  ///
 , cluster(V4_location)

 eststo: logit w4_vacc    $demographics  ///
 , cluster(V4_location)

 eststo: logit w5_vacc    $demographics  ///
 , cluster(V4_location)

 eststo: logit w6_vacc    $demographics  ///
 , cluster(V4_location)

 eststo: logit w10_booster $demographics  ///
 , cluster(V4_location)
 
 esttab     ///   using "SITab2.tex" /// uncomment to save table
, label refcat (notavailable "Availaibility" 1.agecat "\textbf{Age}" 1.incomesource2 "\textbf{Income source}" 1.Q2_demo_marital2 "\textbf{Marital}" 0.Q3_demo_gender "\textbf{Gender}" /// 
0.Q_demo_education2 "\textbf{Education}" 1.Q_demo_religion_enc "\textbf{Religion}" 0.below15 "\textbf{Children}" 0.above65 "\textbf{Elderly}" 0.affectedness "\textbf{Affectedness}" health_satisfaction "\textbf{Health satisfaction}" Q_trust_nat_gvmt_competence "\textbf{Gov. comp.}" 0.Q_attitudes_party_support "\textbf{Party ID}" 0.Q118_infec_infected "\textbf{Covid}" 1.vacc_village "\textbf{Village vacc.}" 1.V6_district_enc "\textbf{District}" , nolabel) ///
indicate("Location fixed effects" = *V4_location ) ///
b(2) se(2) star(+ 0.10 * 0.05 ** 0.01 *** 0.001) stat( N r2_p, fmt(0 2)) replace  nonote
 

 

*************************************************************************************
*END 
*************************************************************************************
